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Abstract: When neglecting capiUarity, two-phase incompressible flow in 
porous media is modelled as a scalar nonlinear hyperbolic conservation law. 
A change in the rock type results in a change of the flux function. Discretizing 
in one-dimensional with a flnite volume method, we investigate two numerical 
fluxes, an extension of the Godunov flux and the upstream mobility flux, the 
latter being widely used in hydrogeology and petroleum engineering. Then, in 
the case of a changing rock type, one can give examples when the upstream 
mobility flux does not give the right answer. 
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Sur le schema aux mobilites amont pour les 
ecoulements diphasiques 
en milieu poreux 

Resume : En negligeant la capillarite, un ecoulement diphasique incompres- 
sible est modelise par une loi de conservation scalaire hyperbolique non-lineaire. 
Un changement dans le type de roche entraine un changcment de la fonction 
flux. En discretisant en dimension un avec une methode de volumes finis nous 
etudions deux flux numeriques, une extension du flux de Godunov et le flux des 
mobiltcs amont, cc dernier ctant largement utilise en hydrologie et en ingcnierie 
pctrolicre. Dans le cas d'un changement de type de roche on pent alors donner 
des cxemplcs ou le flux des mobiltes amont ne donne pas une solution correcte. 

Mots-cles : Ecoulement diphasique en milieu poreux, mobilites amont, lois de 
conservation hyperboliques, condition d'entropie, methode de differences finies, 
methode de volumes finis. 
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1 Introduction 

Under the assumptions that capillary effects are neglected, two-phase flow in 
a porous medium is modeled by a nonlinear hyperbolic equation. In many 
applications, the porous medium is not homogenous. The flow domain has 
to be divided into several subdomains corresponding to different types of rock 
separated by lines or surfaces along which, not only the porosity and the absolute 
permeability of the rock type change but the relative permeabilities also differ. 
This situation is modeled by a single conservation law with a flux function 
discontinuous in the space variable. Numerical methods designed to simulate 
the flow have to be devised to take into account the discontinuities in the flux 
function. 

In this paper, we compare different numerical schemes of the finite difference 
or finite volume type that are used to simulate two-phase fiow in porous media 
with changing rock types. We restrict ourselves to the one dimensional case. In 
the multidimensional case, most numerical methods still use the one dimensional 
flux calculation in the direction normal to the boundaries of the discretization 
cells. Wc also focus on the numerical flux calculation. 

Conservation laws with discontinuous coefficients arise in several other sit- 
uations in Physics and Engineering like in modeling continuous sedimentation 
in clariflcr thickener units used in waste water treatment plants (See |12j , [13j , 
[5]), in traffic flow on highways with changing surface conditions (see [53]) and 
in ion etching used in the semiconductor industry (see [24|). A detailed account 
of the above applications can be found in |26| . Consequently, equations of this 
type have been studied extensively from both a theoritical as well as a numerical 
point of view. 

In particular, for equations governing two-phase flow in porous media, the 
numerical scheme that is commonly used by the petroleum engineers is the up- 
stream mobility flux scheme (sec [3 [551 [50] ) . An alternative finite difference 
(volume) method of the Godunov type based on exact solutions of the Riemann 
problem was presented in Adimurthi, Jaffre and Gowda ([2j). The above paper 
also addressed the problem of prescribing correct entropy conditions at the in- 
terface between rock types and showing the uniqueness of the entropy solution. 
The solutions computed by the Godunov type scheme was shown to converge 
to the entropy solution. The numerical flux developed in [5] is similar to that 
used by Kaasschieter in [21] although written in a more compact form. 

It is natural to ask whether the solution computed by the upstream mobility 
flux scheme also converges to the entropy solution and compare its numerical 
performance with other schemes like the one in [5] and that of Towers ([29j. 
[50]). The goal of this paper is to address these questions. 

In this paper, wc will give an explicit representation of the upstream mobil- 
ity flux scheme for a medium consisting of two rock types and show that the 
numerical flux is monotone. This will help us to obtain estimates in L°°. We 
will then use a suitable modification of the singular mapping technique to show 
that the approximate solution converges to a weak solution of the continuous 
problem. . The key point of this paper is to address whether this weak solution 
is a entropy solution or not. We show by various numerical experiments that 
the solutions computed by the upstream mobility fiux scheme are not consistent 
with the interface entropy condition of [2]. Furthermore, wc construct numerical 
experiments for the case where only the absolute permeability changes. In this 
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case, the solutions given by the scheme do not converge to the entropy solution 
pointwise and differ qualitatively from the entropy solution. The lack of entropy 
consistency leads us to suggest that the upstream mobility flux scheme may not 
be the correct numerical method to simulate two-phase flow in porous media 
with changing rock types and should be replaced by the Godunov type scheme 
develop in 

This paper is organised as follows, In section 2, we describe the equations 
governing two phase flow in porous media with heterogenities. In section 3, we 
summarise the mathematical theory for single conservation laws with discontin- 
uous flux developed in [T], [5] and mention the well posedness results. Section 4 
is devoted to describing flnite difference schemes of the Godunov type as well as 
the upstream mobility flux scheme. The convergence analysis for the upstream 
mobility flux scheme is carried out in section 5. We give explicit representation 
formula for the flux, show that it is monotone and use a variation of the sin- 
gular mapping technique to show convergence to a weak solution. The core of 
this paper is in section 6 where we address the question of entropy consistency. 
First, we consider numerical experiments for the case where only the absolute 
permeability changes and discuss the entropy behaviour of the solutions. Next, 
we construct examples of the case where the relative permeability also changes 
and show that the scheme is not consistent with the interface entropy condition 
of [2]. We derive some conclusions from this paper in section 7. 



2 Two-phase flow equations 

Capillary-free two-phase incompressible flow is modeled by the following scalar 
nonlinear hyperbolic equation 

^ dt ~^ dx 

where (p is the porosity of the rock, 5* = 5*1 is the saturation of phase 1 which 
lies in a bounded interval [0, 1]. The flux function / is the Darcy velocity ^pi of 
phase 1 and has the form 

f = ^1= . ^\ [g + (.91 -52)A2]. 

Here q = fi + (p2 denotes the total Darcy velocity where ipi,£ = 1,2, denotes 
the Darcy velocity of phase i with, for the second phase, 

<P2 = , ^\ [g + (g2 -gi)Ai]. 

Since the flow is assumed to be incompressible, the total Darcy velocity q is 
independent of the space variable x. 

The quantities X(, i — 1,2 may be called effective mobilities. These are 
products of the absolute permeability K by the mobilities kf : 

\e = Kh, e=l,2. 

The absolute permeability K may depend on x and the quantities ki and Xi are 
functions of S which satisfy the following properties : 

ki andAi are increasing functions of S*, fci(O) = Ai(0) — 0, 
^2 andA2 are decreasing functions of 5, ^2(1) = A2(l) = 0. 
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We also shall assume that these functions arc smooth functions of the saturation 
S and so is the flux function /. 

The gravity constants ge, £ ~ 1,2 of the phases are 

dx 

9£ ^ 9P£-r' ^=1,2, 
az 

with g the acceleration due to gravity, pi the density of phase i and z is the 
vertical position of the point of abscissa x. 

Observe that with the above hypothesis, / is a smooth (say Lipschitz) func- 
tion with at most one local maxi with /(O) = and /(I) = q respectively. 

In many practical situations, the porous medium is heterogenous. For ex- 
ample the medium may consist of two rock types separated at the interface 
(x = 0). In this case, the porosity, the absolute permeability and the relative 
permeability change across the interface and the flow is modeled by the following 
equations: 

{H{x)ct>+ + (1 - H{x))r)St + {H{x)f+{S) + (1 - H{x))f-{S)), = 0, 
S{0,x) = Soix) 

(2.1) 

where H is the Heaviside function and the indices - and -|- refer to to the left 
and the right of the interface respectively. The flux functions are given by the 
following. 

f^-T^T±[q+{92~9i)Xtl \f=K^kt. (2.2) 

See Figure 12.11 for notations and the shapes of the flux functions when q < 
0,52 > 5i- Note that the flux functions may also intersect as we will see in the 
numerical experiments of Section [5) 



Rock type I 


Rock type II 







\ \ 


, A , ki , ^2 
< < 1 


cj)+, K+, fc+, k+ 
0< S <1 

f = f+ 


q - 


0- 


9+ \ 





Figure 2.1: Constants, mobilities and fluxes for two rock types 

Equations (|2.ip with flux functions (|2.2p arc a special case for the more 
general single conservation law with flux function discontinuous in the space 
variable considered in [2]. Flux functions /~ and satisfy the following hy- 
pothesis, 

Hi./^,/+ are smooth (say Lipschitz) on [0,1]. 
H2.r(0) = /+(0) = 0, /-(I) = /+(1) = g. 

Hs./^,/^ have exactly one local maximum in [0, 1] with 6^ ~ argmax (/^) 
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and 6^ = argmax {f'^). 

Note that these are precisely the hypotheses for the fluxes assumed in [2]. 

3 The continuous problem 

As presented in the previous section, a change of rock types leads to a single 
conservation law with a flux function discontinuous in the space variable. An 
entropy theory has been developed for equations of the form (j2.ip with the fluxes 
satisfying the hypothesis Hi,H2, H^. We summarize some of the results for the 
benefit of the reader. 

Even in the case where the flux is continuous, it is well known that solutions 
of equations of the above type develop discontinuities in flnite time even for 
smooth initial data. Hence as in the continuous case, weak solutions S of ()2.ip 
are sought and defined as satisfying 

Siptdxdt + j j {H{x)f+{S) + (1 - H{x)f-{S))(p^dxdt 



+ So{x)(l){0,x)dx^O, \fip eC^(RxR+). (3.1) 



It is easy to check that 5 is a weak solution of ()2.ip iff it satisfies in the weak 
sense, 

St + if-iSJh = 0, x<0, t>0, 
St + if+{S))^ = 0, .T>0, t>0, 
S{0,x) = Soix), VxeM, 

and the following interface Rankine Hugoniot condition. 



where 



/+(5+(i)) = riS-{t)) for almost aU t, 



S+(t)= lim S{x,t), S-{t)^ lim S{x,t). 



It is well known that weak solutions for a single conservation law are not 
necessarily unique. Additional admissibility criteria termed as entropy condi- 
tions need to be imposed for uniqueness. For equations of the form (|2.ip . it 
is natural to impose the standard Kruzkhov entropy conditions away from the 
interface a: = 0. These can be stated in terms of the entropy flux pairs which 
are defined as 

Entropy pairs: {(pi,'4'i}i=i;2 is said to be an entropy pair for (j2.ip if ipi is 
convex in [0,1] and 4>[i9) = </);(6')/+'(6l), VaC^*) = 02(f)/"'(f) for & [0,1]. 

Let Sq € L°°(R) be the initial data with < So{x) < l,Va: e K, and let S 
be a weak solution of with < S{x, t) < 1, V(a;, t) e M x R+. 
Interior entropy condition: With 5*0 and S as above, S is said to satisfy an 
interior entropy condition if for any entropy pairs (</3i, '!/'i)i=i,2, S satisfies in the 
sense of distributions 

?^ + 2*(S)<„, vx>0.<>0. 
at ox 

INRIA 



On the upstream mobility scheme 



7 



But for equation (j2.ip , interior entropy conditions like the one above are not 
sufficient to guarantee uniqueness and we need to impose an additional entropy 
condition at the interface. The central issue in the analysis of conservation laws 
with discontinuous flux is the choice of this interface entropy condition. In [5], 
the following interface entropy condition was used. 

Interface entropy condition: With 5*0 and S as above, assume that S^{t) = 
lim S{x,t) and S^{t) = lim S{x,t) exist for almost alH > and define, 

X— »0+ x^0~ 

L = {t>0; S-{t)€{0-,l],S+{t)e[O,9+)}, 

U = {teL; S+{t) = S-{t) = l}U{t&L; S-{t)^S+{t)^0}. 

Then 5* is said to satisfy the interface entropy condition if 

meas {L\U}^0. (3.3) 

This means that the characteristics must connect back to the x-axis on at least 
one side of the jump in the flux i.e., undercompressive waves are not allowed. 
Undcrcompressivc waves i.e., (/^'(<S'+) > 0, f~'{S~) < 0) are unrealistic 
physically as information is not taken from the initial line. 

S is defined to be an entropy solution of (j2.ip if it is a weak solution and it 
satisfies both the interior as well as the interface entropy condition. With this 
concept of entropy solution, the following wellposedness result was obtained in 
[2] , under the following hypotheses on the initial data: 

INi. S'o is such that < 5o(2:) < 1, V xeR, 
IN2. 7V(/-,/+,5o)<C<+oo, 

where N{f~, /+, 5*0) is an estimator of the total variation of the flux function 
evaluated at 5*0. This estimator will be defined precisely below in Section 4. 

We also need the following definition. 
Regular solution. S is said to be a regular solution of (12. ip if the discontinu- 
ities of S form a discrete set of Lipschitz curves. 

The well posedness result is given by 

THEOREM 3.1 Let So satisfy hypotheses IN1JN2 and satisfy hy- 

potheses Hi, H2, H^. Then there exists a weak solution S G L°°(R x M+) of 
\2.1\) satisfying the following, 

(1) For almost all t > and X £ W, S{x^,t), S{x^,t) exist. 

(2) S satisfies the interior entropy condition \3.S\) . 

(3) If S is a regular solution, then S satisfies the interface entropy condition 
\3. 3\) and it is unique. 

Uniqueness is proved by using a Kruzkhov type doubling of variables argu- 
ment. For details, see [2]. Existence was shown by showing that a Godunov 
type finite difference scheme converges to a weak solution and is consistent with 
the entropy conditions. 

We would also like to mention that more recently hypotheses on the fluxes 
have been relaxed considerably to include fluxes of the concave-convex type as 
in |271 H] , and with finitely many extrema as in [S] . Similarly schemes of the 
Enquist-Osher (EO) type have been considered in [Tl[57]. It is to be observed 
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that equations of the type (|2.ip arc special cases of the more general single 
conservation law with discontinuous coefficient of the form, 

ut + {fik{x),u))x =0, u{0,x)=uo{x) 

The wellposedness and numerical methods for this problem are addressed in 
a forthcoming paper [5]. We must also mention that an entropy theory for 
equations of the above type (including a degenerate parabolic term) has been 
developed by Karlsen, Risebro and Towers in [TH] . In [221 EH] , the author de- 
veloped staggered mesh algorithms of the Godunov and Enquist Osher type for 
the multiplicative case i.e. f{k,u) = k{x)f{u). The case with a degenerate 
parabolic term was handled in |18| and the time dependent case in [10] . The 
entropy condition of |19| agrees with that of [2 in many cases but differs in 
certain cases. See Section 6 for a discussion of different entropy conditions for 
equation (|2.ip . For the rest of this paper, we will use the entropy framework 
developed in [2]. 

4 Finite Difference Schemes 

In this section, we present finite difference schemes using for numerical flux 
either the extended Godunov flux pi or the upstream mobility flux used in the 
petroleum industry for the simulation of two-phase flow in porous media. 

Let / be a Lipschitz continuous function. Then the Godunov flux corre- 
sponding to / is given by 

{min f{6) if a < b, 
'^'""frm -f (4-1) 
max f (0) II a > 0, 

The subscript g stands here for Godunov in order to differentiate this flux from 
the upstream mobility flux that we will introduce later. This flux was first 
proposed in [17j and is very popular in the numerical analysis of conservation 
laws. It is based on exact solutions of the Ricmann problem. Let F~ and 
be the Godunov fluxes corresponding to the fluxes and /+ respectively. 

In the case of two-phase flow, the flux functions /^,/^ satisfy hypotheses 
Hi,H2,H3 and the formula (|4.ip can be simplified as follows, 

Pg^o.,b) = min{/"(min(a, 6'_),/"(max(6'_,6)}, 
F+{a,b) = min{/+(min(a,0+),/+(max(0+,6)}. 

These formulas were introduced in [5] and are simpler to implement than the 
general formula (|4.ip . 

Also, following [2], they extend easily to define the interface Godunov flux 
Fg based on the exact solution of the Riemann problem for (|2.ip : 

Fgia,b) = min{/-(min(a,0_),/+(max(^?+,6)} (4.2) 

We remark that the above interface flux coincides with the interface flux ob- 
tained in [21| for which expression (|4.2p represents a compact form, very easy 
to use for computational purposes. It is also easy to check that the interface 
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flux Fg is Lipschitz is both variables, nondccrcasing in the first variable and 
nonincrcasing in the second variable. Note that the interface flux satisfies 

Fg{0,0) = /-(0,0) = /+(0,0) = 0, Fgil, 1) = r(l, 1) = 1) = g, 

but is not consistent. 

Equipped with the definition of the numerical fluxes, we proceed to describe 
the mesh. Let h > and define the space grid points as follows: 

a;_i/2 = a;i/2 = 0, x^+i/a = j /i forj>0, Xj_i/2=jh for j < 0. 

We will also use the midpoints of the intervals: 

'^■'^^^ h for j > 1, Xj = ( ^^7^ ) h for j < -1. 



For time discretization the time step is At > 0, and let t„ = nAt, X = 
For an initial data Sq G we define 

5°+! = - / So{x)dx if J > 0, S°_i = - So{x)dx if J < 0. 

Now we can define the Godunov type finite difference scheme {5*"} induc- 
tively as follows: 



^r+' = S'^-\{F+{S'l,Sli)-Fg{S'l,,S^)), 

S]^+' = S]^-X{Fl{S^,S^^^)-F+{S^_,,S^)) ifj>l, 

S'lV = Sl,-\{Fg{Sl,,S^^)~Fg{Sl2,S^_^)), 

S]^+' = S]-\{F-{S2,S^+:)~F-{S^_^,S^)) if j<-l. 



(4.3) 



Observe that this is the standard Godunov scheme for j ^ ±1, that is, away 
from a; = 0, 

For 5*0 G L°°(]R) and grid length h and At with A = ^ fixed, define the 
piecewise constant function Sh G i°°(M x R+) associated with {5*"} calculated 
by the scheme 



Sh{x,t) = S]' for {x,t) G [xj-i,2,Xj+i,2) X [nAt, {n + I) At), j 7^ 0. (4.4) 

The above Godunov type scheme was analysed in [5]. For this analysis we 
need to introduce 

Nl{rj+,So) = ^ |^^-(5°,^+i)-^^,-(^"-i,^")l 

+ H I ^/ (-^J ' -^i+i ) " ^9 (^3- 1 ' 'S'i ) I 

+ |F,(5°i,5?)-i^-(5°2,^°i)| 
+ |F/(5?,5°)-F,(5^,5?)|, 

A^,(r,/+,5o) -sup7V,?(/-,/+,5o). 

h>0 

It is easy to see that if 5*0 G BV{M.), then Ng{f~J+,So) < C\\So\\bv, where 
C is a constant depending only on the Lipschitz constants of /~ and /+. 

The following convergence theorem was proved. Let M = maxLip{f^ , /+}, 
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THEOREM 4.1 Assume that X, M satisfies the CFL condition XM < 1. Let 
Sq e L°='(M) such that < So{x) < 1 for all x £ R and Ng{f~,f+,So) < oo. 
For h > 0, let X = ^ and Sh be the corresponding calculated solution given by 
(|4.3p . (|4.4p . Then there exists a subsequence hk — > such that Sh^. converges 
almost everywhere to a weak solution S of (j2.ip satisfying the interior entropy 
condition. Suppose the discontinuities of every limit function S of {Sh} form 
a discrete set of Lipschitz curves, then Sh ^ S in Lf^^{R^, Ll^^(M.)) as h — > 0, 
and S satisfies the interface entropy condition. 

The convergence of the scheme was proved by using the singular mapping 
technique which we wiU also use in section 5 albeit with modifications. The 
limit solution obtained was shown to be consistent with the interior as well as 
the interface entropy condition. The key point in the proof of consistency with 
the interface entropy condition was the use of a contradiction argument using s 
test function. The reader is referred to [2] for details. We will use similar ideas 
in the next section. Some numerical experiments involving this Godunov type 
scheme are shown in section 6. 

As mentioned earlier, staggered mesh schemes were proposed in [22], [30] 
and [18] for general single conservation laws with discontinuous flux. In the 
simplified case of a single discontinuity in the flux, the staggered mesh scheme 
of the Godunov type can also be written in the form (j4.3|) by replacing the 
interface Godunov flux Fg with the averaged interface flux Fr{a,b) which is 
the Godunov flux corresponding to the function t = l/2(/~ + /^). This flnite 
difference scheme is analyzed in j29j and is shown to converge for a large class 
of fluxes. Numerical experiments comparing this scheme with a Godunov type 
scheme was reported in |27j . We will also compare this scheme in the numerical 
experiments in Section 6. 

The main objective of this paper is to analyse the upstream mobility flux. 
It is an adhoc flux for two phase flow in porous media, invented by petroleum 
engineers from simple physical considerations, and it corresponds to an approx- 
imate solution of the Riemann problem. The standard upstream mobility flux 
for f~ is given by the following formula: 

F- {a, b) = ^-^l^-. [q + (.91 - 92)Xr] , 

r A7(a) liq + {g, - g,)XT* > i ^ \,2,i ^ I, (4-5) 

\x-{b) if<z+(,g£-.g,)Ar* <0, 1 = 1, 2,^7^ £, 

Similarly, the standard upstream mobility flux corresponding to can be 
defined by the following formula, 

F+{a, b) = ^+^'^^+. [g + (51 - 52)A+1, 

r X+la) liq + {g, - g,)A+* > 0, z = 1, 2, z ^ ^, (4.6) 
A = < £=12 
1 A+(6) if<z+(5,-.g,)A+* <0, i = l,2,z7^£. 



These formulas just say that the mobility Af must be calculated using the value 
of the saturation which is upstream with respect to the flow of the phase (. since 
the sign of the quantity q + [gi — g2)X2 determines the direction of the flow of 
phase 1 and the sign of g + (g2 — gi)X'{ determines that of phase 2. 
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Note that the above formulae are imphcit and have been made exphcit in 
[5]. The flux is shown to be Lipschitz. monotone and consistent in the same 
reference. 

As for the Godunov scheme, an interface upstream mobility scheme needs 
to be defined to take into account the changing rock types. Formulas (|4.5p . (|4.6p 
can be easily extended to obtain the interface flux F{a, b): 

{ A7(a) if q+{gi~ g,)\* > 0, i = l,2,i ^ £, (4.7) 
\* — J £—12 
1 A+(6) ifq+{ge~g,)X* <0, i=l,2,i^l, 

This interface flux preserves the idea of calculating the flux using the phase 
mobilities which are upstream with respect to the flow of the corresponding 
phases. 

Now we define the upstream mobility flux scheme for a medium with chang- 
ing rock types as follows, 

Sl'+^ = SI' - A(F+(S'i", S^) - F(5'!!i, 5i")), 

= S^-X{F+iS^,S^_,,)-F+{S^_,,S^)) ifj>l. 

For So e L°°(K) and grid length h and At with A = ^ fixed, define the 
function Sh G L°°(R x M_|_) associated with {S"} calculated by the scheme (|4.8p : 

5;,(x,t) = 5; for (.T,t) e [a:,_i/2,a;,+i/2) X [nAt, (n+l)Ai), j ^ 0. (4.9) 

We will analyse the scheme (|4.8p in the next section. As for the Godunov 
case we will need a BV type norm which we define as 

N,{r,r,So) = 5] |F-(5°,5^0-F-(5^_i,5°)|+^|F+(50,5°^i)-F+(5,"_i,5°)| 
j<-i j>i 
+ |F(5°i,5?) -F-(502,5°i)| + |F+(5?,^2") -^(5^,5?)!, 

Nirj+,So)^snpNh{rj+,So). 

h>0 

5 Convergence Analysis 

In this section, we show that the solutions defined by (|4.8p . (|4.9p converge to 
a weak solution of ()2.ip along a subsequence as ft. — » 0. We closely follow the 
analysis of [2] and will refer to the above paper for details. We first observe that 
the formulae (j4.5p . (14. 6p . ()4.7p are implicit. The first step is to make them ex- 
plicit. For the interior fiuxes, this has been done in 0. We will give an explicit 
representation of the interface flux. Depending on the ordering of the gravity 
constants, we have to distinguish the following two cases. 

Case 1: gi < g2 

Following [5] , we define the auxiliary quantities for calculating the explicit fluxes 

01 = 9 + (.91 -.92) A2" (a), Si = q + (.91 -.92)A2, 

02 = q + {92 - 9i)Xt{b), 62 = q + {92 - 9i)K- 
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Clearly we have 9i < 62 and di < 82- We have the following lemma for the 
explicit formulae of the fluxes, 

LEMMA 5.1 We can have only the following three cases: 

1. 0<6ii=(5i<6l2 ^ \l = X[{a),X*2= X2{a) 

2. 6I1 (5i < < 6*2 = (52 ^ A* = A+(6),A* = A^(a) 

3. 6ii<6'2=<52<0 ^ At = A+(6),A^ = Aj(6) 

The proof is simple and similar to the case of the interior flux F. Details can 
be found in [55]. This lemma says that just calculating 9i and 82 is sufficient 
to determine the upstream side of the flow of each phase. 
The other case works in the same way. 



Case 2: ,92 < .91 

We define the auxiliary quantities 

di= q + (.91 - .92)Aj(fe), 02 = q + (.92 - gi)Aj'(a) 

and 5i and 82 are as in case 1. Now we have 9i > 62 and Si > 82 and the 
following lemma. 

LEMMA 5.2 We can have only the following three cases 

1. 6li>6'2 = 52>0 Xl=X^{a), A^=A^(a) 

2. 6I1 = (5i > > 6*2 = ^2 ^ A*=A-(a), A* = A+(6) 

3. O>0i^8i>02 ^ A|=A+(6), X^^X+{b) 

Again calculating 0i and 62 gives the direction of the flow of each phase and 
the way to calculate the upstream mobilities. This is easy to implement in a 
code. 

Our goal is to show that the sequence of approximate saturations converges 
to a weak solution of (j2.ip . We begin by stating some of the properties of the 
interface flux. 

LEMMA 5.3 The interface flux F as defined in is Lipschitz in both its 

arguments, non decreasing in the first and nonincreasing in the second argument. 
Furthermore the following also holds 

F(0,0) = /-(O) = /+(0) = 0, F(l, 1) = /-(I) = /+(1) = q. 

Proof. The proof that the flux is Lipschitz is similar to that for the interior 
fluxes in [5] and we omit the details. Also the evaluation of F(0,0) and F{1, 1) 
is easy to check. So let us have a quick pass at the monotonicity properties. We 
consider for example case 1 i.e. gi < .92. 
If < 6*1 = (5i < 6*2 we have 

Ai = A^(a), A^=A2'(a), < (5i = 9 + (,91 - g2)A^(a) < (^2 = 9 + (.92 - ,9i)Ar(a), 
F{a,b) = ^^^^^8i, 9F^^^^) _ {X^na)X2{a)8i + X^{a)iX^y{a)i-82) 



X^ (a) + A2 (a) da^ ' (A^ (a) + X^ (a))2 

Since Aj~,A^ are both positive functions, Aj~ is nondecreasing and X2 is non- 

dF — 
increasing, and < (5i < (^2, we conclude that ——(a, 5) > and F(a,b) is 

oa 
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nondccrcasing with respect to a. Obviously F{a,b) does not depend on b so it 
is nonincreasing with respect to b. 
If 6*1 = (5i < < 02 = we have 

At = A+(fe), A^=A2-(a), <5i=g+(gi-.g2)A2-(a) <0<(52=<?+(,92-.9i)A+(fe), 

pr ^^(^) dF \+{b){\^y{a){-52) dF X^{a){\tnb)5^ 



A+(6) + A^(a) ' da^ ' ' (A+(6) + A^(a))2 ' 95 ' ' ' (A+(6) + A^(a)) 

dF dF — 

Again it is easy to check that -rr—ia, 6) > and -^{a, b) < 0. Therefore F(a, b) 

da . db ^ 

is nondecreasing with respect to a and nonincreasing with respect to b. 

li 9i = 5i < 62 = §2 < Q we have 
AI = A+(fe), A^=A+(6), <5i=q + (,gi-.92)A+(6)<<52 = 9+(ff2-5i)A^(&)<0, 
T( , (A+)'(6)A+(6)5i + A+(6)(A+)'(6)(-52) 



2 



A+(6) + A+(6) ^' db' ' ' {\+{b) + Xtib)y 

dF dF — 

Again it is easy to check that -^{a,b) = and -^{a,b) < and F{a,b) is 

nondecreasing with respect to a and nonincreasing with respect to b. I 

Similar statements for the interior fluxes can be found in [5] . 

Next, we state the CFL condition for stability of the scheme as the following, 

A A/ < 1, 

M = max{ max {^±l/l(Sj, S,+i) l-lIl(S,^i, Sj)}, 

\j\>i.n cla 00 

9F3/2 dF"" dT' dF-^^^^ 

-Q^iS^,S2) - _(5_:,5i) - ^^(5_2,5-i) 

(5.1) 

This type of a condition was explicitly written out in [25| for the case of one 
rock type and a slight modification of it gives the result in our case. Then we 
can prove in a straightforward manner 



LEMMA 5.4 Under the CFL condition i5.1\) . the upstream mobility scheme 
defined by (C73p is monotone. 



We remark that scheme (|4.8p is in conservative form, is monotone, but it 
is not consistent because of the interface flux (some examples are discussed in 
section 6). Hence, the classical theory (see [TTl[m[T3]) does not apply and we 
have to adopt the analysis presented in [2]. The monotonicity of the scheme 
leads to the following discrete contractivity result: 

LEMMA 5.5 Let So,e L°°{R, [0,1]) be the initial data, and let {S^} be the 
corresponding solution calculated by the upstream mobility flux scheme (|4.8p . 
then, 

Y,\S-+'-SJ\<Y,\S-~S-~\ (5.2) 
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Proof: As the scheme (|4.8p is monotone and conservative, this estimate foUows 
by applying the Crandall- Tartar lemma (see [HI). I 

The next step is to obtain estimates in L°° for the approximate solutions. For 
Monotone, Consistent and Conservative schemes , such estimates follow from a 
discrete maximum principle. But for the scheme (|4.8p . the lack of consistency 
implies that the discrete maximim principle is no longer true. Instead, as in [5], 
we can use the consistency of the interface flux at the points and 1 to obtain 
that [0, 1] is an invariant region for the scheme and obtain the following lemma, 



LEMMA 5.6 Let Sq G L°°(M, [0,1]) be the initial data, and let {S^} be the 
corresponding solution calculated by the finite volume scheme (|4.8p . The follow- 
ing holds, 

< 57 < 1 Vj, n. (5.3) 

Proof. Since < 6*0 < 1, hence for all j, < Sj < 1. By induction, assume that 
()5.3p holds for n. Then from Lemma 15.31 we have 

= F_i(0,0,0)<i/_i(^7_i,^;,^]Vi) = ^r' <^-i(l'l'l) = lif ^■<-2> 
= i/i (0,0,0)<i7i(57_i, 57, 5]Vi) = <i7i (1,1,1) = 1 ifj>2, 

= i/_2(0,0,0)<i?_2(5!^2,'5!!i,5r) = 5!!+i<i/_2(l,l,l) = l, 

= H2(0,0,0)<ff2(^"i,5i",5^) = 51'+i<i/2(l,l,l) = l. 

This proves (|5.3p . I 



As pointed out earlier, the key difficulty in the convergence analysis is to 
obtain BV type estimates on the approximations. For a monotone, consistent 
and conservative scheme, such estimates following from the discrete contrac- 
tivity and the translation invariance (see [14j). But in this case, the scheme is 
not consistent and we cannot expect the approximate solutions to be uniformly 
bounded in BV. Rather, the difficulty is circumvented by using the singular 
mapping technique first introduced by Temple in [28j and adapted for schemes 
for single conservation laws by Towers in |29j . The singular mapping was also 
adopted in the convergence proof in [5] . More recently, several modifications of 
the singular mapping have been suggested in [57|, [3] etc. 

The central idea in using the singular mapping technique is to estimate the 
total variation of the approximate solutions under the singular mapping by the 
variation of the fluxes in neighboring cells and use the discrete contractivity. 
This method works well for upwind schemes like Godunov and Enquist Osher 
but it does not work for other types of numerical fluxes like the Lax-Friedricks 
flux. The same is true for the upstream mobility flux and we have to adapt the 
technique to work in this case. We do so by using the idea of chain estimates 
like in j4j . We start by defining the singular mappings. ,we use the following 
standard notation a G M,then a+ = max{a,0},a_ = min{a,0},a = a+ + 
a_, \a\ = — a^. 

The singular mappings are given by, 

Md) - / If-'iOld^, MO) = / \f+\im (5.4) 
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where a G [0, 1] is some number. Note that we use the same form of singular 
mappings as in [2] expect that there are centerered at an arbitrary point as in 
[27] . Now we are in a position to define the transformed schemes for the discrete 
values of the solution. We define them as 

„ ^ f Ms^) if J < -1 ^ r Msr) if J < 1 

\ MS'li) if J>~1 ' ^' \ MS^) if J>1 " 

Like in |27| , we define two sets of transformed variables which enables us to 
simplify the proof to some extent as compared to [2]. Our goal is to estimate 
the variation of the transformed scheme at each time level. For simplicity, let 
us supprcs the subscript n as we are dealing with the same time level. Then 

TVizj) = \zj - Zj+i \ 2^(zj - Zj+i)+ 

In [2] , [57] , this variation was controlled individually in each cell by the flux vari- 
ation across the neighboring cells. For details see lemma (5.4) in [27]. But such 
an estimate relied on the upwind nature of the Godunov flux and is not neces- 
sarily true for the upstream mobility flux as the upstream mobility flux gives 
different answers from the Godunov and Enquist-Osher fluxes when the phases 
are flowing in different directions. Rather, nonlocal variation estimates hold in 
this case as will be explained below. For this, we observe from the definition 
of the singular mapping (|5.4p that {zj — > if and only if Sj > Sj+i. 

Same observation also applies to wj's. We use this ordering of the neighboring 
cell values to define the following, 

Define J_ = {j < —2} and we define some subsets of this set as follows, 

Definition: I_C1 J_ = {i ^ J_: Si < Si+i] is the set of admissble indices. 

Note that this implies that for each i ^ L there exists a unique k{i) such that 

the following holds, 

1. Si < Si^i < . . . < 5',i_fc(i) 

2- Si^k{i)-i < S'i-fc(i) 
We denote the following, 

1. fc(i) = if Si-i < Si and 

2. fc(i) = cx) if Vj < i, Sj > Sj+i 

So there can be atmost one i € I_ such that k{i) = oo. Let be such that 
io = mini. Note that io is not necessarily equal to —2. Now we can define a 

chain as Jj = {j : k{i) < j < i}. With the above definitions, it easy to check 
that J = Uje/J,. 

Similarly denote, J = {j > 1} and we define some subsets of this set as 
follows. 

Definition: I C J = {i <E J : Si < S'i+i} is the set of admissble indices. 

Note that this implies that for each i £ I, there exists a unique k{i) such that 

the following holds, 

1. Si > Si+i > . . . > Si+k(i} 

2- Si^k{i)+i > Si^k{i} 

We denote the following, 

1. k{i) = if Si+i > S, and 

2. k{i) = oo if Vj > i, Sj > Sj+i 

So there can be atmost one i £ I such that fc(i) = oo. We denote the minimum 
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value in / to be Now we can define a chain as J; ~ {j : i < j < k{i)}. With 
the above definitions, it easy to check that J = U^^jJi. 

Equipped the definitions above, we are in a position to state the main lemma 
of this section in the following, 

LEMMA 5.7 Vi G /, the following estimates hold, if i < —2 and k(i) ^ oo 
then, 

+ \F~iS,+,,S,+2)-F-{S„S,+^)\} (5.5) 
for io as defined above, and such that iq ~ we have the following estimate, 

Y^{z,^z,+,)+ < 2M (5.6) 

In case, i happens to be the only index such that k{i) = oo, then 

Y,{zj~zj+i)+ < 2M (5.7) 

Similarly, if i £ I and k{i) =/= —oo, then the following estimate holds 

+ \F+{S,+^,S,+2)-F+{S,,S,+^)\} (5.8) 
for i^ as defined above, and such that i'^ = 1, we have the following estimate, 

E(^j-^j+i)+ < 2Af (5.9) 

And in case i is such that k{i) ~ oo, then we have that 

< 2A/ (5.10) 

Proof. We will only provide proofs for the estimates (|5.5p and (|5.7p . The other 
inequalities follow in the same manner. We have to consider three separate cases 
to check the estimate namely, 
Case 1: < Sk{i) < e~ . 

In this case, it follows from the bounds and the definitions that < S'i < 
. . . < 5fc(j) < 6~ . Hence, one can check that 

E(^.--^.+i)+ = r(^fcw)-r(5.) (5.11) 

From the definition of /, we get that 5*; < S'i+i and >S'fc(i)+i < Therefore, 
using the monotonicity and consistency of the interior upstream mobility flux 
scheme, we get that 

F-{S,,S,+i) < F-{S„S,) = f-{S,) (5.12) 
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and 

F^{Sk(i),Sk{i)+i) > F^{Sk{^),Sk(^)) = f^{Sk{i)) (5.13) 
Therefore by combining the above estimates, we get that 

fiSkit)) ^ f^{Si) = f~iSk{i)) ^ F~{Sk(i)+i,Sk(i)+2\ 
+ F^{Sk{i)+i,Sk{i)+2) - ■ ■ ■ 

Now by using (|5.12p and (|5.13p . wc get that , 

thus proving (|5.5p . Next we consider, 
Case 2: 0' < < 1. 

In this case, it follows from the L°° bounds and the definitions that 6~ < Si < 
• ■ • ^ Sk(i) < 1- Hence, one can check that 

J2{z,-z,+,)+ = .r{s,)-r{Ski^)) (5.14) 

From the definition of /. we get that Si < Si-i and Sk(i)-i < Sk{i)- Therefore, 
using the monotoniticity and consistency of the interior upstream mobility flux 
scheme, we obtain 

F-{S,.i,S,)>F-{S,,S,)^r{S,) (5.15) 

and 

F^{Sk(i)-i, Sk{i)) < F^{Sk(i),Sk{i)) = f^iSk(i))- (5.16) 
Therefore by combining the above estimates, we have 

f~{Sk(i)) - .f^{Si) = .f^{Sk{i)) - F~{Sk(i)+i,Sk(i)+2)\ 
+ F {Sk{i)+i, Sk{i)+2) ~ ■ ■ ■ 
+ ...-F-{S,-uS^) 

+ F-{s,.i,s,)- ns,) 

Now by using ()5.15p and (j5.16p . we get 

Y,izj-zj+i)+ < ^{|F-(5„^,+i)-F-(5,_i,5,)| 

thus proving (j5.5p . 
Case 3:5, < 0- < Sk{i) 

In this case, 3 l{i) £ Jj such that Si < . . . < ^((i) < 0~ < S';(j)_i < . . . < Sk(i)- 
We get that 

^(z,-z,+i)+ = .r{o-)-r{Sk(o) + .r{o-)-r{s^) (5.17) 
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Again by the monotonicity and consistency of the interior fluxes, we have the 
following estimate, 

> F-{e-,e-) = ne-) (5.18) 

Therefore, from the estimates ([STf]) . (|57T8l) . (|5?T2l) and ([5A6l) . we have 

+ F {Sl(i), S';(i)+i) - . . . 

and 

no-) -r (Ski.)) = r(^^-)-i^-(^iw-2-^iw-i) 

+ F (5;(i)_27 — . . . 

+ ■ • • - ^~(<S'fc(i), + 

Combining the above 2 estimates, we get the desired inequality and prove (jS.Sp 
in all the 3 cases. Next, we prove (|5.7p . In case of i being the unique element 
of / such that k{i) = oo. It is easy to from the L°° and Lipschitz bounds that 

i 

— oo 

< M{\9- -i\ + e-) ^ M 

Thus we have the estimate (|5.7p . The other estimates can be proved similarly. I 

We use the above inequalities to show the following variation bound on the 
singular mapping, 

LEMMA 5.8 The transformed sequences are of bounded total variation and 
the following estimate holds, 

max{rF(z;),ry(«;^")} < jN{f-,f+,So) + 2M (5.19) 

Proof: We provide a proof for the sequence {zj^}, the other bound follows in a 
similar way. We 

iel jeJi_ 

(5.20) 
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By adding the chain inequaHties (|5.5p . [531 and l5.7|) . we get that 

Y,{z-~z-^,)+ = 2iJ2 \G{u],u]+,)-G{u]_„u])\ + \G{u"_„u-,~F{u'li,u{)\ 

J<--2 

= ^EK^'-""I+4^'^ (5.21) 

Now by using the discrete contractivity, we get the desired estimate. In a 
similar way, we can get the total variation bound for {w"} and complete the 
proof of the lemma. I 



In order to show convergence of solutions generated by the scheme, we need 
to define the following piecewise constant functions, Let z^,w^ be defined as 
z^[x,t) = z'j,'W^{x,t) = Wj-, V {x,t) e We translate the bounds on the 
discrete values in terms of the above functions in the following lemmas which 
we state without proof (for a proof see P7]). 

LEMMA 5.9 With the functions defined as above and Vt G R+,u;e have, 

max{TF(z'^),T^(u;'^)} < jNh{f~ , f+ ,uo) + AM (5.22) 

LEMMA 5.10 Let Sq e L°°(R, [0,1]) such that N{f-J+,Sa) < oo be initial 
data and let Sh be the co rresponding solutions obtained by the scheme 112.1]) . then 

< Sh{x, t) <1 V (x, t) G M X R+ (5.23) 



J \Sh{x,t)~Sh{x,T)\dx < Nhif~J+,So)i2At+\t-T\) (5.24) 

R 

We are in a position to state our main convergence theorem. The key step was 
to prove the total variation bounds on the singular mappings and the fact that 
the singular mapping is monotone and hence invertible. We have that 

THEOREM 5.1 Assume that the CFL condition is satisfied and the initial 
data So satisfies the hypothesis INi and IN2, Let Sh be approximate solutions 
defined above,then there exists a subsequence (still denoted by h) such that Sh 
converge almost everywhere to a weak solution S of 112.1]) . In fact Sh S in 
^ioc(''^+' -^locC^)) ^'^ ^ ffoes to 0. Furthermore, the limit solution satisfies the 
interior entropy condition i3.2]) . 

Proof.This is the main convergence theorem for the upstream mobility flux 
scheme (|4.8p . The proof follows by the classical arguments of the Lax-Wcndroff 
theorem (See [H]) and the modifications introduced in [2]. We omit the details 
and refer to the above quoted paper for them. 

For any fixed t > 0,we have the BV bounds from the Lemma (|5.9p and by 
using the standard Rellich compactness theorem that upto subsequences (still 
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denoted by /i),we get that t), i(J/i(., t) converge in and for almost all x to 
z{.,t) and w{.,t) respectively. 

For fixed t and for almost all x < 0, then from the convergence we have that 
ip{Sh) z. Note that ip is monotone for x < and we get that Sh{x,t) 
tp~-^{z{x, t)) = S{x, t). Thus for all x < 0,we get that Sh{-, t) converges to t) 
for almost all x < 0. Similarly, for x > 0,we define S{x, t) = il)2^[zh{x, t) and get 
the a.e convergence. Now we use the standard density argument along with the 
time continuity estimate (|5.10p .to get that 

Sh^S eL^MO,T),LUR)) (5.25) 

Once we have the above convergence, we can use the standard arguments of 
the Lax-Wendroff type in order to show that Sh converges to a weak solution of 
(|2.ip .The proof follows exactly as in [27] and we refer the reader to this reference 
for details. Similarly, the consistency of the scheme with the interior entropy 
condition is shown by using the numerical entropy fluxes of Crandall-Majda (see 
[H]). Check 12 for the details. I 



6 Entropy consistency of the scheme 

In the previous sections, we have shown that the upstream mobility flux scheme 
is well defined for two-phase fiow in an heterogenous medium with two rock 
types and that it generates solutions which converge to weak solutions of the 
conservation law (j2.ip and satisfy interior Kruzkhov type entropy condition. 
But, for the solutions of the scheme to be admissible, we have to show that 
they also satisfy the interface entropy condition (|3.3p . As pointed out earlier, 
we remark that this entropy condition is pointwise and essentially amounts to 
the exclusion of undercompressive waves at the interface (a; = 0). In ([5]), it 
was shown that the Godunov type scheme (|4.3p satifies the interface entropy 
condition by means of a contradiction argument. In this section, we investigate 
the question of whether the limit solution generated by the upstream mobility 
fiux scheme also satisfies this interface entropy condition or not. 

It will be shown by means of various numerical experiments that the limit 
solution generated by scheme (j4.8p which uses the upstream mobility flux need 
not satisfy the interface entropy condition and counterexamples are given. In 
other cases that we report, it is far from clear whether the interface entropy 
condition is actually satisfied. In fact, the numerical evidence suggests that 
this condition is not satisfied in the pointwise sense as required in the entropy 
theory of (f^) on account of certain boundary layer phenomena at the interface 
although the entropy condition may be satisfied in a weaker integral sense. 

It is worth mentioning that there are several entropy theories for equations 
of the type (|2.ip like those presented in |2j and in [l9j. It will be shown that 
in some cases, the solutions generated by the scheme (|4.8p satisfy the entropy 
conditions of and in some other cases, the conditions of [5]. 

We will now present the five numerical experiments illustrating five different 
situations. 

Experiment 1 
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In this example, we consider the flux functions given by, 



A+(5) 


= LIS* 




(S) 


= 1.1(1-5") 




= S 


^2 


(S) 


= 1-5 


9i 


= 2 


92 




= 1 





= 1 


q 




= 



As is clear from the above, we are considering that the porosity and the relative 
permeabilities dont change across the rock types and the absolute permeabilities 
only change with /v+ = 1.1 and = 1. The shape of the corresponding fluxes 
/~ and /+ are shown in Fig. \6\ 



f+ : 

0.3 . 




Figure 6.1: Flux functions in experiment 1 



We consider for initial data, Snlx) ~ < ^ .„ „ In this case, since 

^ ^ \ 0.35 if a; > 0. 

the flux functions do not intersect in the interior of the interval (0.1), the entropy 

solution in this case coincides for the entropy theories of [2] and [H] and consists 

of a rarefaction fan joining 0.65 and 0.5 on the left and a steady discontinuity 

at the interface with 0.5 as the left trace and 0.35 as the right trace. Note that 

the entropy solution does not admit underconmprcssive waves at the interface 

as /-'(0.5) = 0. 

We present the solutions obtained by the Godunov type finite difference 
scheme which we term as the Exact Riemann Solver (ERS) and the Upstream 
Mobility scheme which we term as UM. Also, we compute solutions with the 
staggered mesh version of the Godunov scheme developed by Towers in [29j and 
[50] . We term this scheme as AV. The computed solutions are shown at two 
different times and for two different mesh sizes in Fig. 16.21 and Fig. 16.31 
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Figure 6.2: Solutions in experiment 1 with h = 0.1 at times t=0.5 and t=1.5 
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Figure 6.3: Solutions in experiment 1 with h = 0.01 at 1=0.5 and 1=1.5 



In Figure 16.21 we show the numerical results obtained with the mesh size 
h = 0.1 and the CFL constant A = |. As expected at this rather large mesh 
size, the resolution is not high although the ERS is already giving very good 
results with the interface discontinuity being resolved perfectly. On the other 
hand, both UM and AV do not resolve the interface discontinuity well. In 
fact, as seen in Fig. 16. 2[ the left trace of the solution as computed by UM is 
approximately 0.4 which is less than the expected trace of 0.5. This is indicative 
of the development of a boundary layer at the interface a; = 0. Another anomaly 
is the existence of a travelling wave in both UM and AV that is clearly unphysical 
as the solution is constant 0.35 in x > 0. The amplitude of this spurious wave 
is higher for UM than for AV. Both these phenomena indicate that we cannot 
prove that the limit solution generated by UM and AV are consistent with the 
interface entropy condition p.3p in a pointwise sense. 

In order to confirm the above proposition, we reduce the mesh size to 
h ~ 0.01 and show the solution in Fig. 16.31 Again, we see that ERS re- 
solves both the rarefaction and the interface discontinuity very well with little 
numerical diffusion whereas both UM and AV do not match the solution. Even 
with a very small mesh size, the left trace at the interface of the solution com- 
puted with UM is around 0.4 and is well below the required value of 0.5. Also, 
the spurious travelling wave seen before is still present although its magnitude 
has decreased. As stated earlier, the existence of both a boundary layer and a 
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travelling wave forces us to believe that the solution computed with UM is not 
consistent with the interface entropy condition. The same holds true for the 
solutions computed by AV. 

The first numerical experiment that we have presented represents the sim- 
plest type of discontinuity at the interface involving only a change in the absolute 
permeability across the interface. Even in this simple situation, the UM flux 
scheme does not perform as well as ERS and the limit solution obtained by 
it doesnot seem to satisfy the interface entropy condition of [2]. Hence, more 
interesting and complicated behaviour is expected when we consider changes 
in relative permeabilities across the interface. As will be shown in the coming 
numerical experiments, the limit solution computed by UM will converge to the 
entropy solution of [2] in some cases and the entropy solution of [H] in some 
other cases. 

We start with an example where the solution given by UM seems to converge 
to the entropy solution of [2] in the following numerical experiment, 
Experiment 2 

In this experiment, we consider the following flux functions and parameters, 



A+(5) 


= S 




= 2(1-5) 


KiS) 


= 2S 




= 1-5 


91 


= 2 


92 


= 1 





= 1 


q 


= 



In this case, we are changing the relative permeability functions across the 
interface. The flux functions are shown in Fig. El 




Figure 6.4: Flux functions in experiment 2 



Observe that in this case, the flux functions intersect at the point 0.5 in the 
interior of the domain and the point of intersection is undercompressive i.e 
/+'(0.5) > and /"'(0.5) < 0. The initial data are 5o(a;) = 0.5 V a; e M, so 
we start with a state where the light and heavy phases are fully mixed. In this 
case, the entropy solution of [1] is given by the constant state 0.5 connected to 
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the left trace 0.42 by a rarefaction fan on the left and the constant state 0.5 
connected to the right trace 0.58 on the right. Observe that this solution satisfies 
the interface entropy condition jSSl) as /-'(0.42) = and /+'(0.58) = 0. 

As the flux functions satisfy the "crossing" condition of [TH], we can apply 
the Kruzkhov type condition of [12] to get that their entropy solution in this 
case is given by 5" = 0.5. This implies that there is no flow in the medium 
which is unnatural as noticed in [21j . This is one situation where the above 
entropy theories differ and the entropy theory of [2] captures the physically 
relevant solution. We have computed the solutions using all the three schemes 
to obtain the results as shown in Figures 16.51 and 16.61 Fig. 16. 5p shows the 
solutions obtained by schemes ERS, UM and AV with /i = 0.1 and the CFL 
parameter A = 0.125. We show the computed solutions at times t = 1.5 and 
t = 3 respectively. As can be observed in Fig. 16. 5[ the solution obtained by 
AV is the constant state 0.5 in accordance with the entropy theory of [19]. The 
solution computed by ERS converges towards the entropy solution as discussed 
above with a good resolution of the interface discontinuity and some numerical 
diffusion at the rarefactions. The solution obtained by UM shows the same 
qualitative behaviour as that calculated by ERS although the left trace is 0.35 
which is well below the left trace of the solution i.e 0.42. Similarly the right 
trace of the UM solution is 0.65 which is above the required right trace of 0.58. 
This is again indicating the evidence for UM of a numerical boundary layer at 
the interface which was noticed in experiment 1. 

In order to get a better estimate of the boundary layer, we shrink the mesh 
size to h = 0.01 and present the results in Fig. 16.61 
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Figure 6.5: Solutions in experiment 2 with h=0.1 at t=1.5 and t=3 



INRIA 



On the upstream mobility scheme 



25 



- 


ERS 
UM 

\ AV 




\ 
\" 
\ 
1 





■ 


ERS 
UM 

\ AV 




1 





-4 -3 -2 -1 1 2 3 4 -4 -3 -2 -1 1 2 3 4 



Figure 6.6: Solutions in experiment 2 with h = 0.01 at t=1.5 and t=3 



Again, the sohition obtained by AV is the constant state 0.5. As expected 
given the convergence results presented before, the solution obtained by ERS is 
almost exact. Notice that the left and right traces at x = are exactly 0.42 and 
0.58 as in the exact solution of this Ricmann problem showing the high resolu- 
tion of the scheme. The qualitative behaviour of the solution obtained by UM 
is again similar to that of ERS. But the boundary layer at the interface remains 
as the left trace is still below 0.35 (well below 0.42) and the right trace is still 
above 0.65 well above the required right trace of 0.58. This suggests to us that 
the boundary layer remains intact as /i — > and the traces at the interface are 
different from the expected traces, although, the width of this boundary layer 
shrinks with a reduction in the mesh size. This suggests that the limit solu- 
tion obtained by UM converges to the entropy solution of f2] in an integral sense. 



In experiment 2, we considered flux functions where the solutions obtained 
by UM converged to the entropy solution of [2] in an integral sense. The crucial 
point of the previous experiment was that the fluxes intersect in the interior of 
the interval (0,1) and the point of intersection was undercomprcssivc. Next, we 
consider a situation of the similar type where solutions computed by UM seem 
to behave very differently. 

Experiment 3 In this experiment wc consider the flux functions and parame- 
ters given by. 



XtiS) 
91 


= s 

= 2 
= 1 


92 

q 




(1-52) 
1 



Ar(5) 

^2 


= 1.755 
= 0.255 

= 1-52 


+ 0.375 


if 
if 


5 < 0.25 
5 > 0.25 



The flux functions arc schown in Fig. [H) Notice that in this case, and /+ 
intersect at 0.5 and the intersection is undercomprcssivc. This is a situation 
which looks similar to that of the previous numerical experiment. 
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Figure 6.7; Flux functions in experiment 3 



Again wc start with the constant initial data given by So{x) = 0.5 V x £ 

M. 

As in experiment 2, the entropy solution of consists of the constant state 
0.5 connected to the left trace 0.45 by a rarefaction on the left, a steady discon- 
tinuity a.t X = connecting the left trace 0.45 to the right trace 0.54, and the 
constant state 0.54 being connected by another rarefaction to the constant state 
0.5 on the right. As the "crossing condition" is satisfied, the entropy solution 
of [H] is just the constant 5* = 0.5. We show in figure (|6.8p the results obtained 
by all three schemes with h = 0.1 and the CFL A = 0.125. 




-4 -3 -2 -1 1 2 3 4 -4 -3 -2 -1 1 2 3 4 



Figure 6.8: Solutions in experiment 3 with h ~ 0.1 at t=2.5 and t=3.75 

As noticed in Fig. 16.81 the solution computed by ERS approximates the 
entropy solution of [2] while the solution computed by AV is the constant 0.5. 
But what is really surprising is that the solution computed by UM is also the 
constant 0.5. In fact, this example has been constructed in such a way that 
A]" (0.5) A]^(0.5) and A^(0.5) ^ A+(0.5). Hence from the very definition of 
the upstream mobility flux, it is easy to check that the solution computed by 
UM remains the constant 0.5 at all time steps. 
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Thus so far we have shown two experiments involving fluxes with an un- 
dercompressive intersection in which the entropy solutions of [2] and [12] differ. 
In experiment 2, the solutions computed by UM flux seems to converge to the 
entropy solution of [5] in an integral sense, though not pointwise whereas, in ex- 
periment 3, the UM flux gives the constant solution which has been considered 
in literature as unphysical (see [H]) and converges to the entropy solution of 
[H]. Despite similar flux geometry, this inconsistent behaviour of the UM flux 
indicates the difficulties of characterizing the limit solutions computed by the 
scheme. 

The above numerical experiments clearly show that the inconsistent be- 
haviour of the UM flux when the fluxes intersect in the interior of the interval 
(0,1) and the point of intersection is undercompressive. We now investigate an- 
other type of flux geometry in which the flux functions intersect and the point 
of intersection is overcompressive. In this case, the limit solution obtained with 
the UM flux also shows an inconsistent entropy behaviour. 
Experiment 4 In this experiment, we consider the following flux functions and 
parameters, 



A+(5) = 25 


Aj(^) 


= (1- 


S) 


X^{S) = S 




= 2(1 


-S) 


91 = 2 


92 


= 1 




= 1 


q 


= 






Figure 6.9: Flux functions in experiment 4 
.eps 



The flux functions are shown in Fig. [D Observe that / and intersect at 
0.5 and that the intersection is overcompressive i.e /~'(0.5) > and /''''(0.5) < 

0. We consider the initial data So(x) ~ l !p ^ ^ ^ 

^ \ 1/3 if a; > 0. 

The entropy solution of ^2j in this case consists of the constant state 0.66 
connected by a rarefaction to the left trace 0.58 on the left , a steady disconti- 
nuity at the interface between the left trace 0.58 and the right trace 0.42 and 
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the constant state 0.66 connected to the right trace 0.42 on the right. Note that 
the solution is not undercompressive as /~'(0.58) = /+'(0.42) = 0. We remark 
that the above fluxes /~ and /+ do not satisfy the "crossing condition" of [T9] 
and the entropy theory developed in the above reference does not apply to this 
situation. But we can still compute the solutions given by AV as the scheme is 
well defined. We present the solutions in figure (|6.10p . We consider the mesh 
size h = 0.1 and the CFL parameter is A = 0.125. 



0.7 1 , , , , , , , 1 0.7 




°"-4 -3 -2 -1 1 2 3 4 °"-4 -3-2-1 1 2 3 4 



Figure 6.10: Solutions in experiment 4 with h = 0.1 at t=1.5 and t=3 



As shown in Fig. I6.10p . the solution obtained with ERS approximates the 
entropy solution of [5] . Note that the left and right traces are very close to the 
expected values of 0.58 and 0.42. On the other hand, the solution computed 
by both UM and AV is the steady state 2/3 on the left and 1/3 on the right 
which is very different from that of the solution given by ERS. Observe that this 
solution is undercompressive i.e /~'(2/3) < and /■*"'(l/3) > 0. The entropy 
theory of |2] avoids solutions of this type. Also this solution differs from the 
solution of the Riemann problem constructed by Diehl in |12| which in this case 
is identical to the solution computed by ERS. We believe that this undercom- 
pressive solution is unphysical and the right solution is computed by ERS. 

It is easy to show by using that Aj'"(2/3) = Aj^(l/3) and A^(l/3) = A]^(2/3) 
and the explicit definition of UM that the solution computed by UM for all h 
in this case is the steady state with 2/3 on the left and 1/3 on the right. The 
natural question that arises is whether the solutions computed by UM agree with 
that of AV in the case where the flux functions intersect in an overcompressive 
manner. The answer to this question is contained in the next experiment. 
Experiment 5 In this experiment, we consider the following flux functions and 
parameters, 

5052 X+{S) = 5(1-5)2 
1052 X:^{S) = 20(1-5)2 
2 .92 = 1 

1 q =0 

The flux functions are shown in Fig. [51 Notice that in this case, the flux func- 
tions intersect in the interior of the domain at 0.46 and the point of intersection 
is overcompressive. 



A+(5) = 

Ar(5) = 

91 = 
6 



INRIA 



On the upstream mobility scheme 



29 




Figure 6.11: Flux functions in experiment 5 



We consider the following initial data So{x) 



0.8 if .T < 
0.2 if a; > 0. 

In this case, the entropy solution of consists of a rarefaction joining the 
constant state of 0.8 with the left trace of 0.6, followed by a constant state of 
0.6, a steady discontinuity joining the left trace of 0.6 and the right trace 0.32 
and a rarefaction joining the right trace to that of the constant state of 0.2. 
Check that this solution is not undercompressive. The solutions obtained by 
all the three schemes with h = 0.1 and A = 1/32 are shown in Fig. 16.121 The 
solution given by the ERS flux approximates well the entropy solution, even 
with a large mesh size. The solution given by thg AV flux is quite different 
in this case and note that the traces (0.7, 0.22) are undercompressive. On the 
other hand, the solutions obtained by the UM flux are very close to those of the 
ERS flux besides a boundary layer on the right. A further reduction in mesh 
size to h = 0.01 shows that the boundary layer on the right remains and the 
traces given by the AV flux are undercompressive as shown in Fig. ()6.13p . 



ERS: 
UM 





ERS: 
UM 



Figure 6.12: Solutions in experiment 5 with h = 0.1 at t=0.25 and 1^=0.5 
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Figure 6.13: Solutions in experiment 5 with h = 0.01 at t=0.25 and t=0.5 



To sum up about these experiments we observed the following behaviour 
across the interface: 

1. In some experiments (experiments 1,2,5) the upstream mobility flux may 
produce unphysical boundary layers and travelling waves. The traveling 
wave and the width of the boundary layer vanishes when /i — > 0, while the 
heigth of the boundary layer may remain significant. Despite of these nu- 
merical artefacts the solution given by the upstream mobility flux remain 
close to that given by the ERS flux. This suggests that in these experi- 
ments, the solution calculated with the UM flux, even though it does not 
satisfy the pointwise entropy condition ()3.3p . may satisfy some integral 
form of it. For the average flux, depending on the experiment, it behaves 
like the UM flux (experiments 1, 5) or it misses the interface discontinuity 
(experiment 2). 

2. In other experiments (experiment 3, 4) the UM flux as well as the AV flux 
produces unphysical undercompressive solutions (experiment 4) and even 
misses the interface discontinuity (experiment 3). 

7 Conclusion 

In this paper we analyzed the upstream mobility numerical flux for a finite 
difference scheme when a two-phase flow crosses the interface between two rock 
types. This results in a discontinuity in the flux function with respect to the 
space variable. We were able to prove convergence to a weak solution but 
numerical experiments show that it does not satisfy the entropy condition of 

m- 

Most often the solution given by the upstream mobility flux is close to that 
given by the extended Godunov flux but numerical artefacts like boundary layers 
or traveling waves perturb the solution. There are even cases when the upstream 
mobility flux misses the discontinuity at the interface. The solution given by 
the averaged flux is not doing any better. 
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